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1 - INTRODUCTION : OBJET DE L'ETUDE 

Au cours des deux demieres annees, le Groupe d’Analyse Numerique du 
CERT a developpe une methode d’approxlmation de type Galerkin dlscontlnu pour la 
resolution de certains systemes hyperbollques non llnealres dits "symetrlsables" [1,2, 
3). 


Cette methode. Issue de precedents travaux sur les formulations 
variationnelles en entrople, generalise aux systemes conslderes des resultats connus 
concemant les systemes hyperbollques llnealres symetriques (1,3). 

Ces systemes symetrlsables sont tels que le flux de la forme polaire de 
1 entrople, fonctlon des variables entroplques (intensives) admet une decomposition 
en parties corrvexe et concave, globale. 

C'est entre autres le cas des equations d’Euler, pour lesquelles on obtient de 
fagon "naturelle 1 ' cette decomposition, en ecrlvant les equations en variables 
entroplques comme les moments des equations de Boltzmann par rapport a un vecteur 
"de collision" [1, 2, 3), 

Ces proprletes justlflent les autres denominations de la methode de 
Galerkin discontinue comme methode de "Flux Splitting" (decomposition clnetique 
des flux) ou methode de ” 3? - dlagonalisation". 

L approximation "elements finis" de d° O cholsle (equivalant a des volumes 
finis) est Impliclte, et peut etre alsement expllcitee moyennant l'etabllssement d'une 
condition de stability sur le pas de temps, de type "CFL", qui s'lnterprete comme une 
condition de dissipation du schema explicite [4], 

De tres bons resultats, tant theoriques que numeriques, ayant ete obtenus 
en monodimensionnel comme en bidimensionnel, - permettant des comparaisons 
concluantes avec d’autres methodes - [3] - il s'agissait en 1989 d'ameliorer les 
performances de la methode, en precision et temps calcul, par l'elaboratlon, dans le 
cadre bidimensionnel instationnaire, d'une strategic d’adaptation automatique de 
malllage, et eventuellement de domaine, qui soit speclflque a cette approximation 
(exploltant notamment le fait qu'elle peut operer sur des maillages de structure 
absolument quelconque. en espace comme en temps). Au prealable devait se 


t 
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poursuivre l'approfondissement des tests de criteres d'auto-adaptation en 
monodimensionnel. 

D'autre part, l’lmplementation tridimensionnelle de la methode ayant ete 
realisee auparavant dans une configuration simplifiee. en 1989 devalt etre developpee 
totalement 1'approximation 3-D (necessitant l'ecriture des pre-processeurs generaux 
de maillage "automatlque" 3-D) avec test sur une geometrie simple (tout ceci en vue 
d’utlllser ulterieurement l'algorithme specifique d'auto-adaptation de maillage ou de 
domaine. en 3-D). 

Parallelement a ces developpements devalent se poursuivre les travaux de 
recherche fondamentale quant aux approximations de type Galerkin (avec elements 
finis dlscontinus). d'ordre superieur, et a l’etablissement d'une methode de type 
"Petrov-Galerkin" fondee egalement sur la ^-diagonalisation des equations d’Euler. 

2 - INTERET DE L’ETUDE 

2.1 - r^rartdMsHqnp* de la mfithode d e decomposition cln^tlqug dCS flt g 


La methode de Flux Splitting decrite est "naturelle” car d'une part elle 
explolte les proprietes tntrinseques au systeme hyperbollque - dont essentiellement la 
convexlte - d'autre part, elle permet une interpretation statlstique de la viscosite 
numerique (au sens moleculaire) : le systeme peut en effet. par analogic avec 
Boltzmann, s'ecrire comme moyenne par rapport a un vecteur "de collision Jc, de 
1'equatlon hyperbollque scalaire lineaire. 


D'ailleurs dans le cas lineaire symetrique, la decomposition de la fonction 
"primitive" des flux, en partie convexe et concavel^- ^jequivaut 
strictement a la decomposition "classique" des flux selon le signe des valeurs propres 
dejacobien : _ ^ 

J ffwr 


J 


fr c 


<0 


De plus : 

cette methode permet de generaliser en non lineaire la notion de "totale 
diagonalisation M . 



s 


En outre la decomposition en c~ , c~ se fait par un calcul al g CbrldUS 
sans determination de valeurs propres, sans evaluation de max ou min (facilite de 
mise en oeuvre, temps calcul ... ). 

^approximation du type "volumes finis" effectuee sur le schema permet de 
retrouver le schema decentre du ler ordre "classicjue (schema de Courant) sur un 
maillage structure. 


2. 2. a - : 


Le systeme Initial s'ecrit sous forme conservative, au sens des 


distributions: 

111 4 . + 4 ,' qj («n = ° f'A?( t /’ t ^ ' R 

q - ® C<trj X Si. U C £ IR.^ 

Lorsqu’on effectue le changement de variable "entropique" 

*—> tj) = ^S(*i ^ ul est bl J ectif ’ 

S etant une entropie admissible strictement convexe pour (1), on obtient l’ecriture 
equivalente sous forme "intensive" en variable entropique : (non conservative) 

(2) 4 ?*(<?) + M c ly) - o 

On verifie les proprietes sulvantes : 

jVf) Z % 

ou ^f w )]est la transformee de Legendre, ou polaire, de 5(\v). C'est 

une fonction strictement convexe de . 

i - — 

f" (tf) - ** 

ou 5 u * est la polaire, ou transformee de Legendre, du i e flux d’entropie s l 

NB : pour Euler S'" (\jJ ) — ^ S (w ) 

^ composante de la vitesse. 

r- * i -** 

On a alors des relations de polarite suivantes. suffisant a deiinlr ° et o , i' - 4 ^ : 
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( 3 ) 


f 5(wi«fi) + 

S ; (w(Cf)) + ^ Cf. x-.-i. 




On definit alors la fonction 


( < 3 ; = f'T'f ) ~) 


2:*(Cf,S*)= S't'f’W + S 1- *^) /n/.' 


u_ ^ — 

etant generalement la normale espace temps unltalre 
sortante a une frontiere ou llgne de discontinuity interieure au domaine. Alors grace a 
(3) on a : 


% I* = vJ('f) /y U- -H-f'I'f ) 


c'est la somme des flux de la variable non conservative par rapport a 
2.2.b - R61e de la fonction ^ ^ ^ 


On utilise la convention habituelle : 


^ % s*<r ) 


valeurs "interieures 






valeurs "exterieures" 


relativement a V\, 


L’inegalite suivante : 


(4) 


(Sitt-Stf 4 ))* h- ( S(w')-S^)W 

_ S (w e ) . ( { w e - w a ). *W + (^(^) - />U ) 4 


o 


est desormais reconnue pour permettre une description aisee des conditions aux 
limites : sur une discontinuity interieure au domaine (les relations de Rankine et 
Hugoniot etant verifiees) elle equivaut a l'inequation d'entropie. 

Or cette inegalite est toujours strictement equivalente a la suivante : 

151 i-(Cfisr) * 0 

qui n'est autre qu’une inegalite de convexite locale pour ^ en "f .D’ou 1'importan- 
ce de la convexite de £.* dont on montre egalement (cf. [1]) qu’elle est etroitement liee a 
la notion de caracteristiques. 




2.2.C - ; 


Definition : (cf. (3D 

On dira que (2) - ou sa forme equivalente (1) - est "totalement dlagonalisable si 
et seulement si 11 exlste un vecteur &(*],*) dans tR'*' , parametre en fonctlon de (►), E 
£ K fft* . variables d'etat independantes de t et x, tel que la resolution de (2) 

soit equivalente a celle de : f . , , . \ 

/ &<«*.). (Iv +3,(^9)) _ o 


soit equivalente a celle de : f . , , . \ 

/ -Efyi). _ o 

avec la notation : U - X f/ E - 

Alors : si est une entropie strlctement convexe pour l'equation scalaire 

Q 9^ t Q L^) = M^Lpjest le i e flux d'entropie associe et on obtient : 

t ^ ) derivee de la polalre d'entropie, qul equivaut 4 : 

tH » L ' W) 

Dans le cas des equations d’Euler mono, bi ou tri-dimensionnelles < 3 


m-n+2 et : 

I A 

\H 


& 

^ X^X^ 
^X^ X 3 
(^e^p)X 


p*(-v 

^X 2 X 3 

(j*e-*-p) ^ 


& 

| ^XiX 3 

^X Z X 3 

P 4fx; 


p est la pression 

les composantes de la vltesse 
e l’energie totale massique 




est 1'entropie '’physique" 


t 


Alors (par analogic avec les equations de Boltzmann) si on choisit 





ou : 


A ^ ^ - zi- 
5 i 0^0 

lot des gaz polytropiques de 
constante X 


- La fonctlon donne l'equation d’etat (en lamodiflant on peut etendre ces 

resultats aux gaz reactifs). 

Le systeme des equations d’Euler apparait comme systeme des moments 
relativement a X de l'equation scalaire Uneaire. 


On 


Vj 6 K = £ -ro 3 

- (u,, Ui y Ui) - vitesses particulaires - M- 0 




choisit 1^)= C qui permettra de trouver pour entropie { b +■ a ^ ) 
On obtient les expressions suivantes : 

Cas general 


w = 

iTe = 7 # (ip) = v’ s*(ip) 

W = 

1 


— ► , ktp . 1 

k ( n , u ) e dul 

„n J 

- 

£ 


n>o 

" 

/? 


f (<P ) * 


i = 1 , . . . n 


S (tp) = 


i* -*■ 

S (ip) = 


i= 1 , . . . n 


u'iTa = v” s 1 * (ip ) 

ip 


L ( K'P ) 


u L (k.ip! 


f L (1p] 


i» 1 .... n 


i-* K.f. , 
u . k e Ju dn 


S (<P ) = 


^ n>o 

l 

R n 

. n 



r 


. 

kip . 
e du 

n 

J n>o 

* 

R 


dn 


“ 


r 


1 t 




S (») = 



u e d u 


n>o 

. 


i=1 , . . . n 





an 



1 3 



Toutes les quantltes relatives au systeme non lineaire apparaissent 
comme moyennes des quantites correspondantes pour l'equation scalaire lineaire. 



1 1* 


L’interet de cette ecriture reside en ce qu'on diagonalise totalement le 
systeme sans avoir a determiner de valeurs propres. D'autre part, il n'est pas nnn plus 
necessaire d'evaluer, arete par arete, le max ou le min de ( 0 t ^ 


necessaire d'evaluer, arete par 
En effet, on a: 






OU 


_ J f/Ht +jU>0 qui peut s'fecrire : 

o \ 

0) J 


Fxj ^ TUr ^ 


On peut montrer sur l'exemple des equations d’Euler que ces integrales se 
calculent "algebriquement" par un changement de variable ramenant 1 integration a 
un demi-espace : 

supposons n=2 et plagons-nous sur les faces verticales correspondant (en 

pratique nous utillsons des elements droits en temps), 
alors 


r tf ,*) 


” Aj >0 R 


Si on effectue un changement de repere pour se placer en repere normal (U s'agit d une 
rotation). I 

vT (• jjlMx -v-m v 

qj - — 'IT /Vl>e. ■+ 



et un changement de variable sur ^ : 





r 


H 

Lo 

' Ht 1 

1 \ 




'Sl- ^ 


mod 


/ 




J 


>0 





alors : 



*U) -> \ 


(V>o) 


de fagon evldente, 

et ces integrates s’explicitent ,, facilement ,, . Pour obtenir les flux ou derive l’expression 
obtenue. 


Sur chaque element, avec le choix fait pour 
fonction residu ’’erf' ci calculer : 


L* . 


on a une seule 


erf ( x ) 


. -erf I -x 


. f nr 

ou e = — 

0 

Le schema numerique utilise est le suivant : 
n-l : OJ » - XC* est un ele 


est un element flni espace 
- temps interieur a Q = l o,r J X-&- 



1 

t , ] 


1 

n* 1 


1 



1 

n 

d 

1 

-\ 

e 

/ 

\ 


/ 

\ 


/ 

/ 

\ 



\ 


\ 

\ 

t -* 


fonction test et ift fonction approximee. sont continues par morceaux. 
discontinuites localisees sur les 


v u verifiant Su.pp ^ ^ 



i 0 




£ (4 pfqo ,4 ftp) l?2* 

J 0<jjfc V ^ . 

w N r /O i * _ iV", . 


Equations 


c.r. ^ ■*■ ^ L ^ ^ 


choisissant VI* fonction caracteristique de ^ ^ entre autxes on a : 

V2”(<f e V)- V2*'(^V) = o 

*£ u ST 4 - u QSt x j>£d ) 

* r r fdrpq verticales 


faces horlzontales 


/i= 


Alors : 


Soit, si £1 = aire fe) , on a le schema implicite suivant : 

+ At: ( f c ( y )• ^ ) - ( f ( ' ' 

C.I. - faces horlzontales faces verticales 

ou est la solution au precedent pas de temps 

Ce schema peut etre explicite en posant : sur 

zif 4 ) = £ if) -|£ r ( 

a jfep 


ce 


schema explicite peut aussi etre considere comme l ere etape d'unpoint fixe en w, la 
condition de CFL s'identiflant a la condition de convergence. Sur le plan theorique les 
arguments permettant d'apprehender la consistance avec les equations faibles et 
l'inequatlon d'entropie proviennent de la decomposition convexe-concave des flux 

(ref. [11. [4]). 

. Les conditions aux limites utilisees si touche SQ 

(faces verticales) 

On ajoute les termes suivants au second membre : 

* Condition d' obstacle : 

/ 0 



qui equivaut a la nullite des flux des variables conservatives. 

* Condition imposee (a 'Tamont", caracteristiques rentrantes) . 

/ -► <*>77 Uta A 


it L 

valeur imposee ^ = ^6. y 


etant la valeur imposee 


* II n’y a rien sur la frontiere 


"libre" (a "l'aval". caracteristiques sortantes). 


O Q . T '^rithme d’auto- ndap tatlQn du maiUa ^s 


2.3.1 


SClfiflu £ He raffi n gmml SL derafTinement du m affla 


un crltere de ralllnement "naturel" sur une art.* specifique a la methode 
car issu des propr*t6s de l approxin.at.on, de type entropique e. 
est fouml par .estimation de la stabilite du schema, en ImpUclte on en expllclte 
estimation permet d'ailleurs egalement d'etudier l'unic.te d'une solution en 
association avec le choix de la condition limite paroi). 

Aveclechoixdu schema expliciteau temps k on effectue le produit 

scalaire par ^ ( vf interieur "initial" en temps precedent), soit. sur un 

element intdrieur au domaine : , \, 

,1* .r'l- P r At ^ ^- ,z 

On developpe (cf. [Annexe 21) en utlllsant les diverses proprte.es de la 
decomposiUon clnetlque des Dux (proprletes de Z . Z. -2. ) e. on obtient sur 

element ne touchant pas la paroi la relation sulvante : ^ - 


Leterme E represente une approximation decenlree (du fait de 
explication : en implici.e on a .approximation cen.ree.de . inequation d en.ropie : 
c'est done un terme negattf. Le "resldu entropique". K . do,, done necessa.rem 


etre positif ou nul pour que le schema respecte cette Inequation d'entropie, ce qul 
indult ici une condition de stabllite, "de type CFL", sur le pas de temps (majoration 
locale sur Telement, de /\ b ). 

Notons que sur le schema implicite . la stabllite est inconditionnelle, car 
dans le "residu" R on volt alors apparaltre une expression positive du fait de la 
concavite de SL* (inegalite caracteristique de la concavite). 

La posivite de R traduit une dissipation, R doit etre le plus proche 
possible de O (le cas limite R - 0 correspond £ un schema non dissipatif : Tinegalite 
d’entropie devenant alors une equation). 

D'ou le choix de cette expression R caracterisant la dissipation 
comme critere "unilateral" (associe a une normale, sortante a l’element considere ) 
sur une arete. 

2.3.2 - Strategic d’auto-adaptation du malllage et/ou d u domaing 

On choisit : 

-> pour "deraffiner' 1 le malllage, d'associer ce critere ’'bilateralement" aux 
aretes en sommant les contributions provenant des 2 elements adjacents a cette arete 
(en ID ou 2D) correspondant a 2 normales opposees. 

Derafdner consiste done a supprimer les aretes sur lesquelles ce critere bilateral est 
"petit" (apres definition d’une echelle rendant la valeur du critere relative par rapport 
a la solution). 

-> le raffinement porte, lui, sur les elements : il s’agit de subdiviser (on a le 
choix de la strategic de subdivision) tout element sur lequel le critere obtenu par 
sommation des criteres unilateraux sur les aretes est "grand". 

Cette strategic est coherente puisque rimportance de la dissipation 
caracterise la variation de la solution. D'autre part, elle est totalement 
"instationnaire" et permet une evolution automatique du maillage en fonction des 
discontinuites de la solution ("suivi" de chocs ... ) meme dans le cas de phenomenes 
rapides. 



is 


En ce qui conceme les aretes appartenant a la frontiere paroi en espace, le 
critere "unilateral" est legdrement different mais issu lui aussi du schema et de 
l’estimation de stability. Le terme supplemental^, associe a la condition limite paroi, 
dans le schema explicite est le suivant : , . 

-4 . p ft) - d(r ** 2 ' D ' 

paro> ' \o r] 

0 ( 0 ^. paroi etant 1'ensemble des aretes de l'element sur l'obstacle (cette condition 
equtvaut a la nullite des flux des variables conservatives sur la paroi). 

Par consequent, le critere "unilateral" de raffinage entropique prend sur 
les aretes paroi la forme particuliere suivante : j » 


—At 


rme particuliere suivante : ^ . 

p. ** _ ^ 

LO. 0 / 


^ pcu-o* l 0 


-At f (p ^ 
J fa* 


aspects : 


JS Stf^. 

De plus, cette strategic d’auto-adaptation peut etre envisagee sous deux 
-> 1’evolution du maillage dans des frontieres exterieures "fixes" 


-> revolution du domaine et du maillage ; dans ce cas, les frontieres sont 
egalement modiflables (ajout - ou suppression - d'elements - d’arfctes - touchant la 
frontiere exterieure (libre ou avec condition imposee). 

Cette deuxieme demarche convient mieux a priori aux phenomenes 
instationnaires rapides (qui sinon requierent un domaine de calcul important). 


3 - DEROULEMENT DE L'ETUDE 


3.1 - Etude monodimensionnelle 


Le code monodimensionnel integrant l’adaptation automatique du 
maillage a et£ developpe avec le deroulement suivant : 
a un instant donne tn on dipose 

_ d’une configuration du maillage' 


i 
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_ d'une solution ^ sur ce maillage 
le pas de temps est Atw 
on effectue la sequence : 

* calcul des criteres bilateraux sur les "aretes" (qui sont des noeuds) C* 

* derafflnage : etant la borne cholsie pour le critere sur une arete, si ^ A 

alors on supprime 1'arSte A. 

* calcul des criteres sur les nouveaux elements ainsi constitues . &£ 


» raffinage : Eg est la borne choisie pour les criteres sur les elements et si 

C > & alors on subdivise 1' element 
£ ' ^ 

• calcul du nouveau pas de temps AU-M (respectant la condiUon de stabilite) 

* nouvelle resolution par la methode de decomposition cinetlque des flux. 


- Dans le cas d'un maillage frontiere variable", la phase "d'extension au bord 
eventuelle par rajout d'un element s'insere entre le derafflnage et le calcul des criteres 
sur les elements. 


- Apres derafflnage des aretes, on choisit pour tf initial surges nouveaux 

elements constitues. la moyenne, ponderee par les longueurs, des LP sur les 


anciens elements les constituant : 



4* f K 





noeuds supprimes 

- Pour le raffinage des elements, on s’est contente de les subdiviser en 2 sous 


elements de longueur egale (bien des strategiesseraient envisageables a ce niveau). On 
affecte a chacun des 2 nouveaux elements le ^ de l'eiement raffine. On fixe 
egalement une taille minlmale d’element en dessous de laquelle on ne rafflne plus. 


D’autre part, l’ensemble de ces operations a ete effectue systematiquement a 
chaque pas de temps. La encore la strategic est modulable (determiner l’intervafle 
optimal de n pas de temps entre 2 adaptations de maillage). 
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Cependant la complexite du probleme nous a conduits 4 conserver ces 
hypotheses slmplificatrrces pour efTectuer les tests de critere. 

2 - hliU mensionnelle 

II s'agissait d’lmplanter une generalisation bldlmenslonnelle de ce code, 
utlllsant des elements espace de structure polygonale quelconque (proprlete 
essentlelle de la methode, jusque la inexploitee). 

Pour cela on a du realiser un code permettant. de la fa^on la plus generate, 
de mailler des domaines de calcul pour la resolution par elements finis de systemes 
d'EDP du ler ordre llneaires ou non llneaires (cas des equations d Euler). 

Dans le cas de la mecanique des fluides, on utilise done des elements 
destructures qui sont des polygones convexes deflnis k l’etat Initial comme des 
triangles, et qui sont amenes a se transformer suivant 1’utlllsatlon des rafflnements 
et ddrafflnements. On obtlendra alnsl un malllage tres dense au volslnage des chocs, 
permettant de gagner en precision, et au contralre les elements seront tres grands dans 
les zones ou rien naffecte la solution. La flnalite est d'automatlser le deplacement des 
points, le rajout et la suppression des elements (tout ceci blen sur en fonction du 
critere entroplque propre a la methode). Qui plus est. afin de ne pas restreindre k un 
domalne borne le calcul. on choisit de pouvolr etendre le domaine espace-temps de 
fa?on automatique. Toutes ces techniques doivent setendre ensuite au 3D. 


3.2. 1 - Construction du malllage 

^elaboration du malllage necessite la construction des objets d’une part, 
l’ecriture d'outlls de malllage d’autre part. Pour cela un modeleur graphique et un 
mailleur ont ete realises. Les objets 2D sont stockes sous forme d'ensembles de 
domaines plans convexes. bomes par des ensembles de frontleres. (A ces domaines 
peuvent etre associes des codes, servant au calcul suivant le type de probleme resolu. 
Cette structure confere au malllage homogeneite et regularite). 


Chaque frontlere est determinee par l'ensemble des points la constituant. 
qui sont foumis soit par flehier, soit par la "sourls". ou par des crlteres geometrlques 
simples (droites, cercles. ellipses ... ). Les domaines etant numerotes. on affecte a 


22 


. - j-o domalnes contigus. Un domaine est connu par la 

chaque frontiere les 2 numeros des domalnes co ug 

donnCe de . ensemble de ses frontiCres. dans le sens trigonometrique. 

U mallleur s'appuie. au s.ade initial. sue les points 

. . . a, rfAnere les tableaux de connectivtte et la 

frontieres (et cree des triangles) . il g 

numerotation des triangles, aretes et sommets crees. 

Une technique de barycentrage pent etre utilise si le maillage requiert une 
eegutoite : elle consiste a depfccer un somme. au barycentre de 1 ensemble d« noeu 
2 ,'entourent. en s'intenhsant de toucher aux sommets des frontieres codecs. 

Les planches 1 a 4 represented 2 objets et les maillages g^nferts sv 
objets. La concentration des points initiaux dCOntt la finesse du m 
des triangles. 

us Planches 5 et 6 montmnt les ba^centmges effech.es sur ces rnamages. 

3 2.2 - ppfllneme n t anto-adantattf 2 R 

U code de ca.cn. 2D resolvant les Rations M. par * 

Galerkin discontinue. precedemment developpC. traitait le ca 
exclusivement formes d'eiements utangulaires. 

Ls theorte permet.an. .'utilisation d'eiements po.ygonaux -.uelconque. 

_ -f - fait, dans ce sens (on se llmite a des polygones convexes). Ces 
V adaptation du code a ete fai . dans les techniques 

elements polygonaux ont l’avantage d’etre plus soup es 

de rafTinement et deraffinement. 

us etapes de raffinemen. et deraffinement du maillage utlhsent toujoum 

revaluation du reste de 1'inequahon d'entropie du ^“^^“la 
de la variation du flux d'entropie pennet de savolr si on inilue 

solution. 

pour le deraffinement i si la valeur du crl.ere es. ''foible" Unferieure a une 

borne fixee) sur l'arete. on supprime cette arete et 2 , e,e m ent 

l’avoisinaient. Alin de conserver une structure coheren.e on exige 

polygonal qu'll soit convexe. 
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Pour le raffinement : on evalue le critere elementaire (comme indique 
Annexe 2) et si sa valeur est "forte'’ (superieure a une autre borne predefinie) alors on 
rafline 1’element. Le raffinement consiste a appliquer le mailleur sur l'element s il 
n'est pas triangulalre. S'il Test, on le divise en 4 nouveaux triangles par creation de 
points au milieu des aretes. 

Raffinement et deraffinement sont contrbles par le choix de la borne 
correspondante sur le critere : cette borne est fonction de la valeur moyenne prise par 
les criteres sur Tensemble du maillage. 

D'autre part, avec les rafTinements successifs la condition de CFL limitant 
la taille du pas de temps devient tres contraignante. Par suite, ont ete implementes 
une semi-implicitation puis une implicitation locales du schema, afin de supprtmer 
cette restriction rendant les calculs tres couteux (cf. Annexe 4). 

Les figures 7 et 8 montrent des maillages “polygonaux" ou figurent triangles 
et quadrilateres (apres deraffinage). (Theoriquement on n’a meme pas besoin d’un 
"raccordement" au niveau des sommets : les sommets n’ont pas besoin de coincider). 

La. figure 9 represente un objet deraffine entierement puis rafline sur 2 
couronnes interieures uniquement. 

3.3 * Etude tridimensionnelle 


On etudie actuellement la possibility de construire un mailleur 
tetraedrique automatique ou semi-automatique tri-dimensionnel autour d’objets de 
faible complexity, par la methode du front. 

Pour cela est mis en place un mailleur 3D surfacique qui definit alors le 
corps comme un ensemble de faces triangulaires. 

De plus quelques outils de construction graphique 3D sont maintenant 


utilisables. 
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Une etape intermediate a consiste a developper un mailleur 
tridimensionnel utilisable sur des corps axisymetriques. le resultat etant un 
ensemble de prismes et de tetraedres. 

Parallelement. a ete Implantee la version 3D de l'approxlmation Galerkin 
discontinue. 

4 - RESULTATS NUMERIQUES 


4.1 . I JR tube A CHOCS DE SOD MON OniMENSIONNEL 

n s'aglt d’un tube de longueur 1 ferme aux 2 extremites pour une parol : 

1 i 


et initialise par 2 etats (separes par une membrane) 


gauche 



drolte 


0 .AZS 

^.\= o. 

0.2s 


Divers tests ont ete efiectues : 

-> sur un domaine evolutlf : partant de 4 elements (2 de chaque cdte du choc 
initial) de longueur cbtr 5..40"*et effectuant des extensions au bord. (On suppose 
alors le tube "ouvert" aux extremites. on ne peut etudler les reflexions sur la parol). 

-> sur un domaine fixe de longueur 1. correspondant initialement 200 
elements de longueur = 5.»40 


Divers ajustements ont ete necessaires : 

cholx du 4^ initial sur l'element rajoute au bord dans le ler cas ( 
de l’element adjacent ou ^ "exterieur”). 

Un inconvenient apparait a ce niveau : si l’on accepte de "derafflner" le 
bord. on peut passer son temps a aj outer puls enlever au bord le m&me element : par 
centre, si systematiquement on rajoute un element de meme longueur que le volsin. 
on croit exponentiellement 
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- dans le 2° cas, pour le tube de SOD (tres "stifF') la generation d elements 
tres grands est gSnante. On a affecte un Indlcateur aux aretes creees par raffinage. qul 
seules pourront etre deraffinees : alnsl. part ant d'un maillage "pas trop fin , on ne sera 
jamais molns fin. D'alUeurs dans ce cas vlolemment lnstatlonnalre, ll apparalt 
comme necessaire d’avoir Inltlalement une certaine finesse de maillage. 

- Dlverses "interferences” entre raffinage et derafflnage ont lnstaure la 
necessity d'un "retard" au derafflnage alnsl que d'une "avance” au raffinage. 

En effet le calcul du critere expllclte & partlr des resultats precedents provoqualt le 
derafflnage d'elements nouvellement ratlines ! ... D'ou la proliferation du raffinage 
dans les memes zones au temps suivant. Diverses solutions ont ete testees : 


-> modification du critere elementaire de raffinage (cf. [A2D 


aulleude R rr oA A +■ At. 

jL 

on calcule R — Vk. |A| i- -Xjl ■ Cc 

At ^ 


ou. A 4 0 
et CX>/0 V» 


-> "mlxage" entre le critere expllclte et le critere impllclte : 


Part ant du schema expllclte monodimensionnel : 


. si on effectue le prodult par v f > * on obtient la somme : 

inequation d'entrople decentree (expllclte) + critere expllclte = 0 
. si par contre on effectue le prodult par on obtient : ^ 

Inequation d'entrople impllclte + critere en ^ et ^ = 0 


-> simple "decalage" du derafflnage d'un pas de temps (retard impose). 


Les bomes test du critere pour le raffinage et le derafflnage ont du 
egalement etre ajustees pour evlter de provoquer des "escaliers” dans la solution - 
continuite 


Les planches 10 et 11 : donnent un apercu des resultats satisfalsants de 


cette batterie d'essals. 


4.2 - 1 t££ tests menslotmela 


4.2. i . j s. d- * on wdimenstennsl 

u technique dauto-adaptatlon du tnaillage a «t teatee dans le cas du tube a 
Chocs de SOD implements en 2 D lie cas etant en fait monodimensionnei). 

u calcul etant violenunen. instationnalre. on effectue toujours ies 2 Ctapes 

I“": “eTu " toct~lArpenTan, que ie 
deraffinement intervient sur les zones "d’etat constant . 

on obtient ainai une so.ut.on satKfaisante avec une bonne p^ion au 

n^u du choc conune de U detente ,voh couches 12 et 13. oh on a ^nte touted 

in rfr>i - -+ir- H'mi matUaee Initial assez grassier. 


iro iwirc nrises 


4.2.2 - la rfmr^ mclinee 


Un autre cas test instatlonnalre a ete realise : 

U s'agit d'un ecoulement confine sur une rampe inclinee (diddre) 
environ 


d'un angle de 


11° 



L’ ecoulement presente un 


choc fort detache instationnaire. 
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Les figures 14 4 18 represented le maillage au cours de Sexploitation. une 
focallsation sur ,a rone rafflnee. des Iso-valeurs de Mach, pression e. entropie 

massique. 

u visualisation des resultats utilise la nornie graphlque Phlgs. 

^approximation Galerkln discontinue etant constante par element on 
represente les courbes Iso-valeurs pour les resultats e„ colorian. les elements apres 
affectation d une couleur * chaque plage de valeurs. U nombre 
valeurs minimale et maxlmale est llmite par le nombre maximum e co 
Von desire. Cecl permet une Interpretation dlrecte des resultats -bruts sans 
mtervenir aucun lissage nl approximation supplementalre. 


4 2.3 - Qhjpt 2D de forme a y lsvmetriane incidence de 2- 


Un autre calcul a «e effectue en 2D sur le profll symetrique sulvant avec 
une Incidence de 2" (11 s'aglt d’un objet axisymetrlque). 



A 


1 : frontiere libre, "aval” 

2 : frontiere en flux Impose, "amont 1 2 3 

3 : frontiere parol 
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Les figures 19 a 22 montrent une visualisation du malllage et des iso- 
valeurs correspondants. 

4.3 - Rfesultats 3P 

La premiere version simplifiee du mallleur 3D a permis de realiser le 
calcul sur le tube a chocs de SOD sur un cylindre plein. 

Un calcul est en cours sur l’objet precedent mals en 3D avec une Incidence 

de 0.1°. 



La visualisation des resultats sous forme d'iso-valeurs 3D et de coupes 
d'lsovaleurs a ete programmee. On peut grace a un plan que l'on promene sur 1 objet 

visualiser les coupes correspondantes. 

Les figures 23 et 24 representent un debut de malllage 3-D sur une sphere. 

Les figures 25 et 26 montrent l’objet 3D. 
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5 - CONCLUSIONS ET PERSPECTIVES 

— ^ 

La classe des methodes numeriques ouverte par le concept de k 
diagonallsation cinetique des equations hyperboliques, dans le cadre de la mecanique 
des fluides, ne privilegie apparemment pas la methode de Galerkln discontinu 
appliquee avec une approximation d’ordre O. On peut en effet augmenter l'ordre 
d’ approximation. ou fonder une methode de Petrov Galerkln sur le m6me principe (cf. 
[2). [41). Ces methodes ont un point commun : elles conduisent toutes k un contrble de 
l'entropie sur des maillages de structure quelconque. Nous avons effectue des essais 
monodimensionnels en elevant le degre de Tapproximation Petrov Galerkln. Si dans 
le cadre d’un schema implicite les resultats sont corrects (comme le laisse prevoir la 
theorie), il semble difficile de conduire une exploitation conseivant les qualites 
mathematiques de la methode. et des CFL "assez grand" (cf. [5). 

Les approximations de type Petrov Galerkln, quant a elles (cf. [4]) conduisent aprds un 
"adimensionnement" convenable a des resultats de bonne qualite (cf. (6)). 

De toute fa^on, la trop grande dissipation sur les discontinuites de contact, inherente 
a l'ordre peu eleve de l'approximation peut fctre palliee par une adaptation correcte du 
maillage, comme le montre cette etude. 

Cependant, le caractere explicite de la methode et la dimension des 
elements spatiaux necessaire a une bonne approximation conduisent a des couts de 
calcul eleves (pas de temps trop petits). Cest pourquoi, dans ce cadre, il est necessaire 
de passer & une formulation "seml-implicite" qul conserve la mdme precision spatiale 
en etant plus performante (stabilite inconditionnelle). 

Nous avons pourtant rencontre une difficulte d’ordre numerique (en cours 
de resolution) : le systeme local non lineaire a resoudre dans le cas de la semi- 
implicitation, est mal conditionne quand le Mach tend vers rinffni (en pratique Mo£ 
V). En effet, il provient de la stationnarisation d’une fonction convexe, mais dont l'une 
des valeurs propres tend vers 0 quand le Mach tend vers -h». Lorsque ce probleme sera 
resolu, soit par un changement d’entropie, soit par un conditionnement local du 
systeme elementaire, il semble que rimplementation en 3D condu.’se a des resultats 
interessants tant en precision, qu'en robustesse et en performance. 

D’autre part, cette methode est immediatement generalisable (cf. [4], [5]) 
aux melanges hyperboliques de gaz reactifs, tels qu’on peut en rencontrer en 
hypersonique. 11 semble de plus, que l'interpretation numerique de la diagonallsation 
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via les equations de Boltzmann, sur chaque espece peut egalement permettre de 
prendre en compte plusieuis equations d’energie. c'est-a-dire de trailer un cas de 
melanges reactifs de gaz hois equilibre thermodynamique. en conservant la structure, 
les methodes de resoluUon et les justifications mathematiques des equations d'Euler 

des gaz parfaits. 
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ANNEXE 1 


EXPLICITATION DE LA DECOMPOSITION CINETIQUE 
DES FLUX EN DIMENSION 1 
(pour la dimension 2. cf [3)) 


*n= 1 dimension d’espace 

/ A \ ( Ax, vitesse 




1 


j parametre reel positlf 


A +■ - _d 


'2, ' ' s i nZT 

vecteur "de col!ision M 

solt ‘cl - = £ 

• la polalre de l'entropie 5(w> = -£- est: 






s *(<^ (-^0 ^ 

et on a 3| 

o r- 


I{&) - j &~z~ <dvj 

e *1 d*) 




c (<ri - _i\ 

r U 1 


c(t)= ^ =^J 0 ^ dv ) : 

X (2.) ~ vl CU) = >fi f c'^ dv) - f]f 

•* r\ 


done 



-*.<i 

M 

« 

M 

‘ 

V 

^-t 

^t. 

r^l 

M 

3l 

l 

i * 

w 







I! 



CO 

cn 
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ANNEXE 2 


ETABLISSEMENT DE L'EXPRESSION DU CRITERE ENTROPIQUE 
D ADAPTATION DE MAILLAGE 
LIEN AVEC LA STABILITE. 


* On ralsonne en bidimensionnel, sur des elements triangulaires en espace, cylin- 
driques en temps (prismes espace (temps). 



^ est constante sur 


element espace temps interleur au 
domalne 

hi jta-H f , constante sur l'element 




Le schema s'ecrlt (sans Integration par parties, sur un element interleur) 

V7^r * -/ rtJ -*> 


{ fV) = o ■ 


>*. — ► j 

Si on multiplle par on obtient 




1 ( vz YvV ) - vz ^ = ° 


* En dissociant temps et espace, tenant compte de l'approxlmation et en explicitant le 
schema sur les faces "laterales", on obtient (et c'est une estimation de la stability du 
schema) 




terme temps 


terme ne portant que sur l'espace 
(n normale espace ici) 


: aire de 





At : pas de temps - tvi 

• On utilise les propriety sutvantes (cf 111, [3]) ^ 

S't'l')- w*^). Y - 5 s ^ ^ ^ 

" 2 ~ 0?. 7" J , tv X) 

r J 

fl ^ ; ' • j — \ ■ 

_.J- • • « -n vf*-7Z*(CJ‘-n) T- Z'T^^'f -2 Cf/ m 

-:un r 75 -\ 4 M J « 1 + ^^'- 


v-\ / ' - — , 

£1 %■:?.*) - v **■<*■• - r ®* ) r !?:' 1 

» ' — \ S‘Tr?: Q ' ^ ^ (<J) Z? \ 


on a done pour* schfcna -J** f ( ^ ) _ Zffa * ) ) 

‘ ^ 0 


On a a one puui ic 

0* v Ab 


+ At.f 

J &1 ■* — •" " 

terme spatial 


VZ. ) 


I — c 


Le terme temporel s' 

4 (« 


. « (r ^ f 1 7>- 


_ ^ +«h f“ v 


i)-s-a i -))A^- s( <‘\ 


R — o 


Flnalement on voit apparaitre : 


3 S 


<r s <S ^ j s (^fs) 


1 ~ 

r <3*. A + Ab 2_JL . (U 


ou 

et 


«k' _ 

R _ <5- f S*(^)- s*(^ - (Y^-^VS*^; )J . 

R - L J it j (riijV) - 1 ft 4 # +ir-ri v 2 )) 

l’integrale sur Ml s’exprtmant 
comine so i nine sur les aretes de Sc 

>*— i 

Xv est la longueur de l'arete i _*J -*C\ \ 

est constante sur l'arftte i de . 

On reconnalt enEun terme negatif ou nul car : E < 0 est l’lnequation d'entropie 
decentree du fait de l’explicitation du schema. 

Le residu R doit etre posltlf. pour qu'il y alt stabillte (cela correspond a une dissipa- 
tion) et le plus proche possible de 0. 

or A 4 0 par convexite de 3 . 

et Vv : parconcavitede 2"' (inegalites caractdrlsU^) 

done si Vb (U -O alors necessairement A— ° et R = ^ • 

R est choisi comme critere de raffinage sur un element. (Sur une arete, le critere 

"unilateral" est la valeur de C ) 

i? _ err A •+• Atr Gj ^ o » z. 

or A : 

car ce terme est du ler ordre en At 

c i_ \ c* GT/-* fc) _ Alr^ (*, ) . VS 0 ^ ^ 

Jen effet 5.^1*,^)= * . Y l ' 1 i AP 

_ 

« ( Xj h-i) = 4 1 K t») _ a(t ) tO(|) ) 

<3t 

La positivite de R est done equivalente a la condition de stabillte locale sur 
l’element (de type CFL) suivante : 

y ^ 

- At f( -t- At Z-tCc ^ 0 ou At ^ 2-^ G- 


^iri 


K 


i-r 



b 0 


En impliclte on aurait le schema : >. 

Jk \ i~D K V7* <^*(10* ^ LP» <r • ( Cii 'l ♦. S n ) 


vJ(^) _ 7 * S*(^ ).^f* _ s* ^ ^ * ) 


iA ' 

1 H-l 


0,t( S’lfVS^f) - w(f).^ 

«. At j (zr^-zi^,«))- 4 

<£c 3^ 

o- G£ ( s*t ^ ) + SC^“ ) -» (Cf 4 Cf ^ ) VS" (<f -S*(^)-S(^f ) 

i- At f - At( 

\ J Jk 


Finalement :E + R = 0 _ . 

00 e- <*(s^VS(^)) + At 


«cf 


a_Ctf S €m~P 


£ 40 eJ’l’inequatlon d'entropie €S^' 

K. <*(s*$~)-S*$~) + (T -1 ^) VS'^) 

+ At [ 2*i^- )-ytf i l ») -(? ) 

•'S.T a ^ — — ' -- rt 


J Sz< 


Qfi>\ 


icucc seuJJ. . 


et dans ce cas on pent ecrire globalement (espace/temps) 

t L - 2' (T, 0 

M. ' 

►<A -♦? 


4 


Z'lif >1-1 «V 




NE : Dans l’un ou l’autre cas, la stability globale est etudiee en sommant sur les 
elements Tequation E + R = 0 
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ANNEXE 3 


expression des criteres de raffinage 
et de deraffinage en monodimensionnel 


Le terme residuel de l'lnequation d'entropie constitue un critere "unilateral" sur 
me artte (= un point en monodimensionnel) du maillage : 

. — _ V . _ — k J v y — * dl 


ae arete (= un point en monodimensionnel) du mauiage . — v \ ^ n 

r (zTi?V)^V) ^WV)- I 

::'T ^ 

D6rafflner eonalste a aupprtnrer une artte aur laquelle un tetere est -petit" ; pour cek, 
on calcule un critere Mattel par artte. en aonanant les contrtbuUons unllatenr- 
les" venant des 2 elements adjacents. 

• 


m n etant fixee on obtient : 

etque Z * , * ) - (qf, *0 ♦ 2* Wj* ) 

= ]^ z*(^v)- Z\q to + ( (vz'^^)rVZ* ( ^ ) 


ichant que 
et que 


. > o 

monodimensionnel : 


^ : d-1 ^ je 

J,, i-i t** 




F ) = V Z ('f , ^ 

Fp Ff)- H’,*) 


_ F"(^) — F"( ’“f*"') + ( +fmf^ 


si 


i- 4 , 


gi - '"{dunfe 


4 4 




* Pour rafflner un element : on le subdivise (en 2 lei dans un premier temps) si le critere 
R calcule par sommation des contributions "unilaterales" de ses aretes, additionne 

de la condition initiate est "grand . 

Dans ce cas : ((ft Cft*) V£ * 

c^cJU, -Z (3*v).,Z ) V2 - 1L ? ' 

/.. [ z* c ^' 

C ""' ’J«feW)- j r i«. f'ifemeJ'fe- 

= Jio. A J^^uecu- ote. 

v- ( X' (u> ) = -Z r*- } -o. — \ 

* i^-c^/c). - vs * ( ^ =- wlvf) - 

Soit : 

R =i AtT(Cgauche + Cdroite) + Cinitial fflkrA F* (Cj^N 

et Cdroite * FI I 3J" ) - F~ iV') 

cecl pour i < k < me. -j fe _, 

sik= 1 on modifie Cgauche ou _JL -J* # j -u 

si k = ne on modifie Cdroite ou ^ ^ ° 
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ANNEXE 4 


SEMI-IMPLICITATION LOCALE DU SCHEMA 


Le schema implicite, inconditionnement stable, s'ecrit : 



Avec rintegration par parties, equivalant & rajouter le terme nul : 
- on a Tecriture equivalente : 


sur un element 
interieur au 
domaine espace- 
temps 




^approximation etant constante par element espace et tranche de temps, on detaille 


de la fagon suivante : 

<T t f W 


i juivaiuv, . j * w w | 

(Cf 4 *)- 5T (<)) At 


fr\ , 


ou 


est l’aire de. Tg. 


on eflectue une seml-explicitation (ou semi-implicitation, par rapport au schema 

— * ^ 

explicite conditionnellement stable) sur le terme en , d'ou : 

w tf) +. Ah j_ 72'^ V) = w 7Z’ 


cela correspond a une relaxation de type Jacobi. 


Ce schema est encore conditionnement stable mais de fagon beaucoup moins 
restrictive que l'explicite. II s’agit d'un systeme non lineaire impliquant tous les 


voisins. 

On peut cependant lineariser ce systeme par pseudo-Newton ; en efTet on peut l’ecrire : 


P l^) = VG F« 

^ aj-wf#*) -At ( 72*'l^y) 
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et 


&(<f J ) = OfcS’f? a )+- At 


G est une fonctlon ^onvexe strictement monotone, ce qui assure que le systeme 
VG ( ^ ^ admet une solution unique et peut etre resolu par une 

methode iterative de pseudo-Newton : 


initialisation 


X. = % 

x- = x~- Jr(xn.F(x-) 

etsl f*= [ tc u 






£ 1 +. comporte le hessien de 




On calcule les derivees par differences divisees : \ 

^ J, frt.fr) -fit*,- 

n- " a 

Cette methode est utilisee localement sur chaque element, ce qui est assez couteux 
malgre le gain sur le pas de temps, bien qu'on se limite au calcul d’un seul Jacobien par 
element (qui est reutilise au cours des iterations). Cette semi-implicitation perm 

dejade multiplier le At par 10. Elle accentue un peu la dissipation. 

On pourrait egalement resoudre par gradients conjugues ; U s'agit de minimiser 

^ de la fa ? on suivante : 

VG(X / ") 4 - ° . par exemple 

minimise chaque ligne : 

p _ b u- 

‘ *X£ • HakG-lX*) it" 

dependant, pour se liberer completement de la contrainte de CFL sur A . 

•galement resolu le schema implicite par une methode de point fixe : 


u'= VS (X") 


on 


4 / 

gui converge vers la soluUon du schema ImpUclte dnitlalement on prend If/ ). 

A chaque etape du point fixe, on resoud par le pseudo-Newton, et cecl sur c aque 

element espace du maillage. 

- A t ne dolt cependant pas tendre vets — stnon .e terme de dlsstpat.on tendratt 

Tten qu^ pas de temps utdise putsse an etre trts grand (dans le cas de condgumOot* 
asymptottquement statlonnalres du moms,. Cette methods es, malgre tout couteuse en 
temps CPU (& cause des 2 methodes ltdrattves imbriquees). 

Jusau'a _ 3 . eUe donne de bons resultats, bien que lente. 

^ , - M ~ ~ £ \ in methode de Newton se met a 

Mais a des Machs plus eleves (on a essaye H*_S)lamet 

_ j au * mtervlent dans le Newton est 

diverger. En fait, la matrice hessienne de 4- q 

tres mal conditionnee dans ce cas (or elle depend du )• 

D faudrait alors effectual un changement dentropie pour que cette matrice ne depe 

plus du Mach. 







Figure 


2 : Maillage 1 


Fiaure 3 : Definition du maillage 






Figure 4 : Mai 11 age 






Figure 5 : Barycentrage du raaillage 






Figure 6 : Barycentrage du mailiage 






Figure 7 : Deraffinage du maillage 1 . 



Figure 8 : Deraff inage cu raaillag 



Figure 9 : Raffinage sur une partie du maillace 







Figure 10 


Tube a chocs de SOD ID 


Maillage auto-adaptatif 


Initialement 50 elements de longueur 


T = 0.1797 


entropie massigue C ^ 


densite 


x 10 


vitesse 


pression 







Tube a chocs de SOD ID 
Maillage auto-adaptatif 

initialement 200 elements de longueur 


T = 0.1920 


entropie massique 


densite 


x 10 


vitesse 


pression 








Tube a chocs de SOD 


<<2-D» en cours de raff inage. 












.787 



Rampe inclinee M» 
lignes iso-pression 




Figure 20 : Corps de rentree 2D M<*> = 8. 
incidence 2° 

Maillage adapte a t = 4.2 













•saaqds aun.p jnoqne a£ squanraTP.p ..sqonoo,, aun.p uoxqonjqsuoo : 








